General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



>-3 to 


J ASA-CE- ! 7026 5) A PHTHCD FCF EIAGKOSING 
FACE EAEAMETEFS G E I N G GECSTPTICKPEY 
ELLITE IMAGFFY A N E A EC C K C A F 1- I A Y E R ECEEI 
s » Thesis (Pennsylvania State Ur.iv*) 

■ 0 p HC A 0 7/M F AO! CSCI 04E G3/47 


F 83-24050 

Onclas 
! 1506 



The Pennsylvania State University 
The Graduate School 
Department of Meteorology 


A Method for Diagnosing Surface Parameters 
Using Geostationary Satellite 
Imagery and a Boundary-Layer Model 

A Thesis in 
■.Meteorology 

Arthur C. Polansky 



Submitted in Partial Fulfillment 
of the Requirements 
for the Degree of 



November, 1982 




The Pennsylvania State University 
The Graduate School 
Department of Meteorology 


A Method for Diagnosing Surface Parameters 
Using Geostationary Satellite 
Imagery and a Boundary-Layer Model 


A Thesis in 
Meteorology 
by 

Arthur C. Polansky 

Submitted in Partial Fulfillment 
of the Requirements 
for the Degree of 

Master of Science 
November, 1982 


I grant The Pennsylvania State University the nonexclusive 
right to use this work for the University's own purposes and 
to make single copies of the work available to the public on a 
not-for-profit basis if copies are not otherwise available. 


Arthur C. Polansky 


11 


We approve the theala of Arthur C. Polanaky. 


Date of Signature: 

SyrfJ? fttz- 

4 Oct /frl~ 


DU- tf gib 


Toby Nr Carlson, Profeaaor 
Mete 


of 


Jotyh 

M 




A. Dutton, Professor cf 
teorology, Chairman of the 
Department of Meteorology 



id K. Blackadar, Professor of 
Meteorology 


iii 


ABSTRACT 

A method for diagnosing surface parameters on a regional scale via 
geosynchronous satellite imagery is presented. The method is a modifi- 
cation of one described by Carlson and Boland (1978) . Moisture 
availability, thermal inertia, atmospheric heat flux, and total 
evaporation are determined from three infrared images obtained from 
the Geostationary Operational Environmental Satellite (GOES) . The 
output fields compare favorably with fields generated by the 
Carlson et al. (1981) method, which uses HCMM data. 

Three GOES images - early morning, midafternoon, and night - 
are obtained from computer tape. Two temperature-difference images 
are then created. The boundary-layer model of Carlson et al. (1981) 
is run, and its output is inverted via cubic regression equations. 

The satellite imagery is efficiently converted into output-variable 
fields. All computations are executed on a PDP 11/34 minicomputer. 
Output fields can be produced within one hour of the availability of 
aligned satellite subimages of a target area. 
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1.0 INTRODUCTION 


1.1 THE PROBLEM 

Remote sensing of the terrestrial environment has provided earth 
scientists with a wealth of data, and is beginning to yield economic 
benefits. Noteworthy are the economic and research benefits derived 
from weather satellite imagery and LANDSAT imagery. Weather satel- 
lites are valuable sources of timely meteorological data for opera- 
tional forecasters, and are becoming significant data sources for 
operational numerical weather prediction (NWP). LANDSAT has 
provided valuable data to economic geologists, cartographers and 
lapfi'v'j!?^ planners (Short et al. , 1976). 

The principal operational weather satellites are equipped with 
two types of sensors: visible-light detectors and detectors sensi- 

tive to the thermal emission of the earth (longwave IR) . LANDSAT 
has detectors for visible light and solar infrared reflected from 
the eqrth (shortwave IR) . Weather satellites may thus be able to 
study ground properties invisible to LANDSAT by detecting variations 
in surface infrared emission. 

The longwave IR sensors on the Geostationary Operational 
Environmental Satellite (GOES) and Heat Capacity Mapping Mission 
(HCMM) satellite measure effective surface blackbody temperature. 
Several soil properties can be deduced from surface temperature and 
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temperature change during the daily heating and cooling cycle. A 
simple energy-balance, equation for the earth's surface can be written: 


S(l-A) + If * H + LE + If + G 
o o o 


( 1 - 1 ) 


where : 

S = shortwave flux incident on the surface (insolation) 

A ■ albedo (visible and shortwave IR) 

If ■ longwave IR flux from atmosphere, incident on surface 

If * longwave IR flux emitted by surface 

H q * sensible heat flux from surface to atmosphere 

L » latent heat of evaporation 

E q =» evaporative flux (mass of water) 

G q * sensible heat flux into the ground 
The longwave IR emissivity of the ground is assumed to be 1. 

S is determined by geometry, the solar constant, and atmospheric 
moisture and turbidity. A can be estimated for most land-use cate- 
gories. The satellite measures a flux related to If. If the 
remaining terms can be estimated, soil moisture availability and 
thermal inertia can be determined. 

Soil moisture availability (M) can be defined as the ratio 
of actual evaporation to the evaporation that would occur from a 
flat surface of pure water. M influences evaporation (E q ) , since M 
is a direct parameterization of the amount of water available to 
evaporate. Thermal inertia (P) , is defined by the equation: 


p = (C A ) 
g s 


1/2 


( 1 . 2 ) 
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where: 

C * ground heat capacity (volumetric) 

© 

X * thermal conductivity of the soil 
s 

P is a measure of the ability of the ground to store heat during 
the day and release it at night. G q is determined largely by P. 

In this work, a numerical model is used to provide values of 

II, H , E , and G . M and P are input parameters of the model, and 

o o o 

are varied over a range of possible values. Model output is then 

compared with satellite data, and actual M and P values are diagnosed. 

H and E values are also obtained, at no extra cost, 
o o 

1.2 HISTORY 

The use of the longwave I.R detection capability of weather 
satellites for remote sensing of soil characteristics is in its 
infancy. Wetzel and Atlas (1981) have developed a technique for 
diagnosing soil moisture from morning surface temperature change 
and windspeed. Dejace and Megier (1979), following Rosema et al. 01978), 
have mapped soil moisture, ground thermal inertia, and daily evapora- 
tion using aircraft measurements of longwave IR and HCMM imagery. 

Carlson et al. (1981), following Carlson and Boland (1978), diagnosed 
moisture availability, thermal inertia, and the surface energy 
budget over two cities, using HCMM data. The above approaches are 
similar to one another. A computer model of soil and lower atmos- 
phere is run, using a variety of soil characteristics as initial 
conditions. A set of forecast surface temperatures is generated. 
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Observed data is then processed, and the model is inverted mathe- 
matically, to yield soil characteristics from the observed 
temperatures. Price (1980, 1,982a, 1982b) has adopted a different 
approach. Surface temperature data are incorporated directly into 
the basic energy-budget equation for the surface, yielding estimates 
of soil parameters and evaporation. 

This thesis is an extension of Dodd's work (Dodd, 1979), 
reported by Carlson et al. (1981). The method is similar to the 
Carlson and Boland (1978) method. Dodd used HCMM imagery to obtain 
surface temperatures; my work compares HCMM and GOES imagery. Dodd 
studied a small-scale phenomenon characterized by large differences 
in surface parameters: the urban heat island. I extend Dodd's 

technique to studies of regional-scale patterns of surface moisture 
variability, which I relate to antecedent rainfall differences and 
crop moisture index variability over agricultural regions. 

The method of diagnosing soil parameters from satellite data 
developed by Carlson is described briefly in Section 2. Several 
modifications were made in this procedure when GOES imagery was 
incorporated. These modifications and their implications are dis- 
cussed in Section 3. Of particular importance are the model 
sensitivity studies presented in Section 3.6. Two case studies are 
presented in Section 4. Conclusions and suggestions for future 


research are collected in Section 5. 
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2.0 THE CARLSON /DODD METHOD 

Quoting from the abstract of Dodd (1979): 

A flexible analysis system has been developed which 
combines high-resolution satellite-derived radio- 
metric ground temperature information with output 
from a numerical model of the boundary layer to infer 
the spatial variation of thermal inertia (P) , 
moisture availability (M) , and the surface energy 
budget. 

The system consists of three basic phases: 

1) Data acquisition and preprocessing 

2) Numerical boundary-layer simulation 

3) Production and display of output fields 

Satellite images are received on magnetic tape. Utility 
programs LONLAT and REGGIE are run to select sub .mages of identical 
size, shape, and geographical location from a day-night HCMM image 
pair. Also, the raw satellite intensities are converted into 
effective surface temperatures, using atmospheric moisture values 
obtained from soundings. 

A one-dimensional numerical model developed by Carlson and 
Boland (1978) is used to simulate atmospheric and ground response 
to insolation. The model and its use in Dodd's study is completely 
described in Section 2.1 of Dodd (1979), and will not be repeated 
here. 

Model simulations are run for 16 ordered pairs of (M, P) values. 
Output from the model at simulated satellite overpass times is 
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passed to a regression routine. Regression equations for M, P, H o , 
and total evaporation (E) are produced, using surface temperature 
at satellite overpass times as the predictors. Biquadratic regression 
equations omitting the cross-product term were used in Dodd (1979). 

The regression equations derived from the numerical model are then 
used to convert the satellite subimages into output images of M, P, 

H q , and E. The output images are then plotted or displayed on a 
graphics terminal. For a complete discussion of this method, refer 
to Section 2 of Dodd (1979). Henceforth, this procedure will be 
abbreviated CD. 

Note the following features of CD: 

1) CD utilized HCMM data. Only two image times per day were 

used: approximately 1400 and 0200 local time. 

2) The postprocessing phase produced biquadratic equations 

in midafternoon temperature (TD) and night temperature (TN) . 
Higher-degree regression equations were not used. 

3) CD was Implemented on a mainframe computer, Penn State's 
IBM 370/3033. The numerical model and regression routines 
in particular had to be submitted as batch-processing jobs. 
Slow turnaround and the mechanics of batch job submission 
rendered interactive use of the model impossible, 

These are the principal features of CD which have been modified to 
produce the system to be described in Section 3, the GOES method. 



3.0 THE GOES METHOD 


3.1 THE PROBLEM 

CD uses HCMM polar-orMter data to diagnose M, P, H q , and E 

over small areas. Dodd (1979) studied the urban complex. 

Kocin (1979) estimated M over a small rural watershed. Carlson 

and DiCristofaro (1981) discussed the applicability of diagnosed 

H q in estimating plume spread over urban and rural areas. In 

3 4 2 

these studies, the target area was on the order of 10 to 10 km . 

For other applications, such as regional crop assessment, target 

5 6 2 

areas on the order of 10 to 10 km may need to be studied. 

Either CD should be demonstrated successfully on this larger scale 
(henceforth called ’regional scale'), or another method should be 
devised to handle regional-scale domains. 

As Kocin (1979) noted, studies of the development of M and 
P anomalies with time are difficult using polar-orbiter data. Due 
to its orbital parameters, HCMM provides a day /night image pair 
only once per 16 days for target areas in midlatitudes. The situa- 
tion is similar with TIROS N, but not quite as restrictive. Cloud 
cover at either image time can render the affected image pair 
unusable. Thus, diagnoses of M, P, H q , and E derived from HCMM 
images, using the CD method, can have no better time resolution 
than 16 days. For studies of the time-dependent behavior of soil 
parameters, a different data source is desirable. 
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CD was restricted to two image times per day, since HCMM never 
scanned an area in midlatitudes more than twice in a day. However, 
the HCMM overpass times, approximately 1400 and 0200 local solar 
time, were close to ideal, given the restriction of only two image 
times per day (Carlson and Boland, 1978). These optimum times 
occur near the times of maximum and minimum surface temperature, 
respectively. If observations of surface temperature were available 
more frequently during the diurnal cycle, better diagnoses of M, P, 
H q , and/or E might have been possible. Wetzel and Atlas (1981) 
suggest that the time derivative of temperature may be better cor- 
related with M than the actual temperature. The possibility that 
additional data could improve the CD method is certainly worth 
investigating. 

In this thesis, a modification of CD is presented which uses 
GOES images instead of HCMM images for surface temperature data. • 

The greater spatial and temporal availability of GOES data addresses 
the problems cited above. The rationale behind the use of GOES 
imagery is described in Section 3.2. Theoretical justification is 
presented in Section 3.3, and practical considerations are con- 
sidered in Section 3.4. Section 3.5 discusses the image-differencing 
feature of the GOES method. An extensive series of quality and 
sensitivity tests were run on the GOES method; the results are 
reported in Section 3.6. 
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3.2 WHY USE GOES IMAGERY? 

CD is theoretically valid only when: 

1) HCMM provides a day/night image pair of the target area 

2) nonadvective conditions prevail 

3) skies are clear 

Obtaining all three conditions simultaneously over an area of interest 
can be a major challenge. Obtaining all three simultaneously on a 
regular basis over a period of months is a forlorn hope. Even 
fulfilling condition (1) alone is difficult. Dodd (1979) used a 
night/day image pair in both his urban complex cases. Kocin (1979) 
used an image pair (June 9-10) in which the night overpass time was 
36 hours before the day overpass time — a departure of 48 hours 
from the modeled night overpass time. GOES imagery is available 
every day; HCMM image pairs are available only once per 16 days. 

If GOES data are used, the assumption that nonadvective conditions 
last more than a day is never necessary; if HCMM data must be used, 
this assumption of stationarity must often be made (as in Kocin, 1979). 
If GOES data are used, no one day is crucial. If clouds cover the sky 
or advection is significant on any particular day, the next day's 
data can be used. With HCMM data, clear skies and nonadvective con- 
ditions on a day/night overpass-pair day are crucial: the next 

day /night pair is not available for over two weeks. Even under the 
restriction that we must use images at 1400 and 0200 local standard 
time (LST), GOES imagery provides sixteen times the chance of meeting 


the three criteria listed above. 
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GOES imagery is produced every 30 minutes. This greater 
availability of imagery adds another dimension of flexibility to 
the diagnostic technique. GOES imagery can be used to diagnose 
relevant variables even if: 

1) skies are partly cloudy 

2) the clouds do not cover enough of the sky long enough to 
invalidate the energy budget calculations of the model 

3) the clouds are moving fairly rapidly 

HCMM imagery would almost certainly fail under these conditions. 

With either source of imagery, areas obscured by clouds at any 
image time used must be discarded. With HCMM, only two image times 
exist, and both are needed for CD. Soil parameters could not be 
diagnosed at any point where there were clouds at either HCMM over- 
pass time. With GOES imagery, there is flexibility in choice of 
image times. If a crucial spot is covered by clouds at one image 
time, a previous or subsequent image could be used (but see Sections 3.4 
and 3.6). Also, if hardware or software failure causes one GOES 
image to be lost, a subsequent or previous image could be used. Such 
a failure could cause a HCMM case to be discarded. 

The greater availability of GOES data allows a more fundamental 
change in CD. There is no need to restrict the diagnostic method 

to two images or to any particular image times if GOES data are used. 

48 

There are 48 GOES image times per day; in principle, there are 2 
possible sets of image times which could be used. 
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Incorporating GOES data into the CD diagnostic method provides 
the following advantages: 

1) more days with images at convenient times 

2) greater flexibility in selecting image times. 

3) the option to use more than two images in the method 
These advantages must be weighed against the disadvantages of using 
a geostationary platform instead of a platform in low earth orbit: 

1) poorer resolution (larger pixel size) 

2) greater distortion of the image, due to the curvature 
of the earth 

For regional-scale studies, the advantages outweight the disadvan- 
tages. For high-resolution, small-scale studies, HCMM data 
probably should be used, if available. The greater flexibility 
resulting from the use of GOES data may be important even for small- 
scale studies. 

3.3 THE GOES METHOD: THEORETICAL JUSTIFICATION 

CD used two images per day because only two were available. 

Wetzel and Atlas (1981) implied that the choice of image times 
and particularly the number of images used in CD may not be optimal. 
With the greater data base provided by the GOES satellite, the restric- 
tions on image times encountered in CD are removed. Consideration 
of the effects of M and P on the daily march of surface temperature 
suggests that three or four images per day should produce superior 
diagnoses. The reasoning leading to this conclusion is outlined in 


this section. 
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M is a measure of the amount of water available for evaporation. 
It therefore plays an important role in the partitioning of energy 
between sensible heat flux and latent heat flux during the day. A 
high M value implies more evaporation, and less sensible heat to be 
distributed between ground and air. Carlson and Boland (1978) and 
Wetzel and Atlas (1981) showed that the effect of M on surface 
temperature is greatest during the morning. In the morning, available 
moisture will be evaporating, and the mass of the air beneath the 
decaying nocturnal inversion will be small. A relatively small 
change in the ratio of latent to sensible heat can produce a rela- 
tively large change in the morning temperature rise, since the heat 
capacity of a thin layer of air is small. In the afternoon, the 
boundary layer reaches its greatest thickness. Temperature rise is 
slower, since the thick layer has a larger heat capacity. The 
magnitude of the temperature rise becomes less sensitive to M. At 
night, evaporation ceases or becomes negative; M becomes indeterminate 
and its effect is inconsequential. 

P is a measure of the ability of the ground to store heat 
during the day and release it at night. The effect of P on surface 
temperature should be most pronounced during early evening. As 
insolation decreases and the surface begins to cool, atmospheric 
convection subsides. The atmosphere tends to decouple from the 
surface, particularly if the wind is light. The surface cools more 
rapidly than the atmosphere, since the surface is a better radiator 
than the air. During the evening, stored heat from below must supply 
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most of the energy radiated from the surface to space. A large P 
value implies relatively large amounts of heat will be stored in 
the soil, and that there is sufficiently large thermal conductivity 
to allow the heat to reach the surface from below. Large P Implies 
slow cooling of the surface. Later at night, the downward thermal 
emission from the relatively warm atmosphere partially compensates 
for the radiative loss from the surface. If wind and stability con- 
ditions are right, an increase in wind speed above the surface layer 
will occur. Low Richardson-number conditions will develop, producing 
turbulent transfer of heat to the surface. P remains important in 
determining surface temperature throughout the night, but its rela- 
tive importance decreases as the night progresses. 

The analysis of the effect of P on surface temperature during 
the day is more difficult. Certainly, if energy is released at 
night, it must be stored during the day. As the sun rises and 
insolation increases rapidly, P must have a strong effect on surface 
temperature. However, as surface temperature rises, the atmosphere 
becomes convectively coupled to the ground. This coupling enhances 
evaporation and atmospheric heat flux. Each flux increases at its 
own rate as the morning progresses. In the early morning, when 
surface temperature is most sensitive to variations in the magnitude 
of heat sinks, latent heat flux and atmospheric heat flux are both 
significant and both varying. Ground heat flux tends toward a 
minimum during this period (Sellers, 1965). The effect of P tends 
to be masked by the effects of M and H . CD model simulations 
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(Carlson and Boland, 1978) show that P is more closely correlated 
with night temperature than with midafternoon temperature. 

3.4 GOES METHOD-. PRACTICAL AND NUMERICAL CONSIDERATIONS 

As stated in Section 3.3, M is closely related to morning 
temperature rise; P is probably closely related to evening tempera- 
ture fall. With only two image times available, it is not possible 
to express both these rates of change. The success of the CD method 
under the data-availabiXity restrictions of HCMM is encouraging. 

With GOES Imagery, we have the freedom to use more images and dif- 
ferent combinations of image times, to capture any portions of the 
diurnal march of temperature that might be useful. Three image 
times could be used: 

1) early morning, just aft<;, sunrise 

2) midafternoon 

3) near mi.dnight 

Alternatively, four images could be used to isolate the two tempera- 
ture changes which should be most closely related to M and P: 

1) early morning 

2) late morning 

3) midafternoon 

4) near midnight (late evening) 

In principle, all 48 GOES images for a given day could be used. 

There are several fundamental problems with the use of multiple 
imagery in CD-type models. These problems are linked to uncertainties, 
errors, and approximations characteristic of all numerical models and 
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remote-sensing instrumentation. In this section, we consider the 
practical and numerical restrictions on the use of imagery in the 
diagnostic model. 

Although the HCMM and GOES satellites both have high-quality 
radiometers, both have limitations in resolution and accuracy. 

These limitations are more severe for GOES. Both must view the 
earth through an atmosphere containing a variable amount of water 
vapor and a variable vertical temperature profile. Corrections for 
water-vapor absorption are applied to the raw satellite data, 
greatly reducing errors but not eliminating them. The effect of 
thin cirrus clouds on the quality of satellite imagery is difficult 
to assess. 

The satellites detect thermal infrared radiation (longwave 1R) 
incident from below. A narrow band of wavelengths is measured to 
minimize atmospheric interference. The emission is converted to a 
temperature by inverting the blackbody radiation equation 
(Goody, 1964, Sec. 2.2, Eq. 2.35-2.37): 



(3.1) 


where: 

B-^ = Planck function at wavelength X (energy flux per unit of 
wavelength) 

C x = 2TThc 2 

C 2 = hc/k 


h = Planck’s constant 
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c * speed of light 
k - Boltzmann's constant 
9 -> thermodynamic (absolute) temperature 

This temperature is assumed to be the surface temperature. True 
surface temperature is the temperature of the earth's surface: the 

temperature of the infinitesimal layer of condensed-phase matter 
forming the lower boundary of the atmosphere. True surface tempera- 
ture, unlike the conventional surface air temperature, is difficult 
to measure. A thermometer placed in a standard Instrument shelter 
1.3 meters above the ground measures the surface air temperature, 
by definition. The reading will be fairly independent of the con- 
struction of the thermometer and the shelter, within reasonable 
limits. If a thermometer is placed in direct sunlight on the surface 
of the earth, a temperature can be recorded. However, this tempera- 
ture will be Influenced as much by the properties of the thermometer 
a s by the temperature of the surface which the thermometer obscures 
and possibly deforms. In a vegetated region, surface temperature is 
hard to define; its measurement may be an intractible problem. One 
could debate where or what the 'surface' is: is the surface the top 

of the plants in a particular area, or is it at some level in the 
canopy? Hand-held IR radiometers can be used to measure surface 
temperature. They reveal significant small-scale variability. 

Cooper (1981) reports that a single point on a vegetated surface 
will h$ve a range of 'observed' temperatures: observed IR temperature 

is a function of the orientation of the instrument and configuration 
of the canopy in many cases. 
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The boundary-layer model used in CD Is comparable to many one- 
dimensional models, but is nevertheless only an approximation. It 
is one-dimensional; the atmosphere is three-dimensional. The model 
is discrete in space and time; the atmosphere is continuous. The 
atmosphere is not limited by approximate physics, and is not subject 
to truncation error and closure assumptions. The model output is 
converted into regression equations, to invert the model mathematically. 
These regression equations represent the model well, but not perfectly. 

The problem of diagnosing the surface parameters reduces to 
the problems encountered in using: 

'regression equations instead of actual model runs 
•an imperfect numerical model 
•imperfect satellite measurements 

•’surface temperature' - a quantity which may not be definable, 
and therefore, not measurable 

In addition, we must assume that the 'surface temperature' used in 
the numerical model is the same physical quantity as the 'surface 
temperature' observed by the satellite. The errors and uncertainties 
inherent in these problems have significant numerical implications. 

These uncertainties place several restrictions on the choice of 
image times and analysis techniques which may be used in a GOES-based 
diagnostic model. The principal restrictions involve the Interval of 
time between image times, and the total number of Images which can be 
used . 

To obtain an estimate of the midmorning time-derivative of 
surface temperature, it is necessary to use a pair of image times 
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bracketing midmorning. In theory, the smaller the time interval 
between the two images, the better the approximation to the deriva- 
tive. In practice, with uncertain data and uncertain model output, 
a short time interval between images will produce a very poor estimate 
of the derivative. Refer to Appendix 1, Section 1 for details. Image 
times several hours apart are needed, to guarantee a large temperature 
change between image times. Errors in the estimate of the derivative 
introduced by the nonlinearity of temperature as a function of time 
must be accepted. Figure 3.1 (from Wetael and Atlas (1982) after 
Schmugge et si. (1978)) shows that morning temperature rise and 
evening temperature fall are close enough to linear to permit finite- 
difference approximations to derivatives using data times several 
hours apart. 

The temperatures at the image times become che predictors in a 
multivariate r^olynomial regression equation foe each diagnosed 
variable. Numerical constraints on the number of image times and 
the permissible degree of the regression equations are discussed in 
Appendix I, Section 2. To avoid overfitting the model output, it 
is necessary to restrict the total number of regression coefficients 
and the degree of the regression equations. 

To maximize the actual temperature differences between images 
while minimizing the total number of images used, the following 
image times were selected: 

1) early morning ; 

2) midafternoon (near daily temperature maximum) ; 

3) around or after midnight. 
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Preliminary testing was done, simulating the use of two images (CO) 
and four images, as well as simulating the above image times. The 
three-image simulation was superior to the two-image simulation. 

A four-image simulation, using late-morning and mid-late-afternoon 
instead of the midafternoon time, provided no Improvement over the 
three-image simulation. Since numerical difficulties and cost 
increase as the number of images increases, the three-image method 
was adopted. 

3.5 IMAGE DIFFERENCING 

The three images used as input to the GOES method are dif- 
ferenced, producing two temperature-change fields: 

DTD * (midafternoon temperature) - (early-morning temperature) 
DTN - (midafternoon temperature) - (night temperature) 

There are three reasons for performing this difference, as will be 
explained in this section: 

1) to approximate the derivative of temperature 

2) to reduce numerical errors 

3) to correct for large-scale temperature variability across 
the target area 

It has been shown that M and P are strongly correlated with 

3T 8T 

-gr . DTD is a crude approximation to at midmorning. DTN iq a 

3T 

crude approximation to in the early evening. In Section 3.4, it 

was shown that these crude approximations are adequate and that 

8T 

better approximations may be unattainable. If is to be used to 
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diagnose M and P, itself must be estimated from satellite data. 

St 

Numerically, DTD and DTN are probably the best estimates of :£■=• 

at 

obtainable from satellite data at the present time. 

Using DTD and DTN as predictors, a bicubic regression equation 
of the following form can be used: 

X - a 00 + a ll DTD + a 12 D ™ + a 21 DTD 2 + a 22 DTN 2 + a 31 DTD 3 + a^DTN 3 

(3.2) 

There are seven regression coefficients. Were three temperatures 
to be used, ten coefficients would be needed for a cubic regression 
equation. It is unclear whether the additional coefficients could 
be added to the regression equation without overfitting the data. 

Using differenced fields appears to be the best way of taking 
advantage of a higher-order regression equation. 

When three quantities are mapped into two, some information is 
inevitably lost. In this case, when three temperatures are mapped 
into two differences, the information discarded contain a potential 
source of error. The information discarded when computing DTD and 
DTN is contained in the mean: 1/3 (T^ + T 2 + T^) . 

The target areas for the studies presented in Section 4 are 
several hundred km in longest dimension. Over distances that large, 
it would not be unusual to find systematic temperature differences, 
even in relatively uniform air masses. The numerical model used 
to simulate the boundary layer over these areas is one-dimensional. 
Horizontal temperature differences cannot be modeled. Only one 
initial surface temperature and only one initial sounding can be 
used. If the target area is dominated by a relatively uniform air 
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mass, the vertical structure of the lower atmosphere will be similar 
throughout the region. Soundings from various points in the target 
area would have similar profiles when plotted on skew-T/log-P charts. 

In the absence of advection, any two points in the region with iden- 
tical surface characteristics should experience the same temperature 
rise and fall during a 24-hour period. However, one of these points 
could be a few degrees cooler or warmer than the other point through- 
out the period. Consider an extreme example. Let the surface 
temperature be 300 K at point A in the region, and let the surface 
temperature be 305 K. at point B in the region, several hundred km 
from A. Assume identical surface characteristics at the two points. 
Also assume clear skies and negligible advection throughout the region. 
Assume that the soundings at tha two points have similar shapes in the 
lower 200 mb: sounding A is uniformly 5 K cooler than sounding B. 

Let both points receive the same total insolation. Will the surface 
temperature difference between point A and point B vary significantly 
with time? 

The term in the surface energy budget equation (Eq. 1.1) most 

4 

directly dependent on surface temperature is It - aT . G^ and H q 
depend on vertical temperature gradients in ground and air, respec- 
tively. E q depends on surface temperature, vertical temperature 
gradient, relative humidity, and vertical moisture gradient. We 
will ignore the absolute-temperature dependency of E q here, with 
ample justification. The 5 K temperature difference between points A 
and B represents a difference of (305/300)^, or a factor of 1.068 in 



23 


total radiation emitted from the surface. If we make the reasonable 
assumption that vertical gradients of temperature and moisture are 
the same at A and B, we can ignore effects due to G , H , and E . As 
a rough estimate, assume that DTD is directly proportional to energy 
surplus at the surface. A 7 percent change in It might then produce a 
7 percent change in DTD. If DTD ■ 35 K at point B, then DTD at 
point A would be 1.07(35) * 37.45 K, a difference of about 2.5 K. 

It is important to note that It is not the only term in Eq. 1.1 
which depends on absolute temperature. It depends on an effective 
average temperature of the atmosphere. If we assume that this average 
atmospheric temperature (T ) at B is 5 K higher than at A, It at B will 

cl 

be larger than It at A The compensation will not be perfect, since 
I a < T in most cases, and the effective emisslvity of the atmosphere 
differs from the effective emisslvity of the surface. The compensa- 
tion will be significant. In a numerical simulation, using the 
boundary- layer model, lowering all input temperatures by 5 K lowered 
all output temperatures by about 5 K. DTD and DTN changed by about 
0.5 K. 

Since a temperature gradient on the order of 1 K/100 km produces 
a variation of only about 0.5 K in modeled surface temperature 
change in a regional-scale study, the effect of large-scale tempera- 
ture gradients on the GOES method is unimportant. If there is a 
large enough temperature difference across the target area to produce 
unacceptable errors in modeled DTD and DTN due to radiative effects, 
the nonadvective assumption will be invalid and the GOES method cannot 


be used. 
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3.6 GOES METHOD: SENSITIVITY AND ACCURACY STUDIES 

In any 24-hour period, there art) 17296 distinct ordered triplets 
of GOES times. Even if the choice of GOES-time triplets is restricted 
to a morning image, midafternoon image, and a night image around or 
before local midnight, there are nearly 1000 possible triplets. 

This section describes the selection of GOES time triplets which can 
be used in the GOES method, subject to the data quality and numerical 
restraints considered in Section 3.4. A method for evaluating GOES 
triplets is presented in Part 1. The results of the evaluation are 
tabulated in Table 3.1 and are discussed in Part 2. 

3.6.1 An Evaluation Method for GOES Time Triplets 

The best way to evaluate the performance of a numerical model 

* 

is to test it against actual data: 

'Gather accurate data for several contrasting situations 
•Simulate each situation numerically 
•Compare numerical results with observations 
This procedure could not be followed for the GOES method. M, P, H q , 
and E are difficult to measure. It would be logistically impossible 
to measure them on the regional scale. Even a set of local measure- 
ments in several selected locations would have proven economically 
infeasible. The difficulties in measuring (or even defining) sur- 
face temperature have been mentioned in Section 3.4. The numerical 
difficulties described in Section 3.4 must also be considered. 


A model simulation which provides satisfactory results with high- 
quality ground-based temperature measurements might fail when 
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Table 3.1A 

Sensitivity Testing Summary: M 

Worst 


GOES Times Output Storage 



(LST) 


Value 


3M 

Error 

r 

1 

2 

3 

M-.5 

8DTD 

3DTN 

M-0.5 

Percent 

06 

14 

23 

.5541 

-.1593 

.1029 

.1043 

98.0 

07 

14 

23 

.4726 

-.1387 

.0674 

.0808 

95.3 

08 

14 

23 

.4972 

-.1383 

.0135 

.0595 

90.2 

09 

14 

23 

.4996 

-.1992 

-.0191 

.0861 

94.1 

10 

14 

23 

.4901 

-.2714 

-.0389 

.1216 

92.6 

11 

14 

23 

.4997 

-.3689 

-.0333 

.1603 

89.4 

07 

09 

23 

.6034 

-.2739 

.1428 

.1652 

89.5 

07 

10 

23 

.5839 

-.2269 

.1238 

.1393 

91.4 

07 

11 

23 

.5764 

-.1930 

.1079 

.1195 

93.4 

07 

12 

23 

.5406 

-.1709 

.0938 

.1048 

95.2 

07 

13 

23 

.5510 

-.1520 

.0885 

.0948 

96.1 

07 

14 

23 

.4726 

-.1387 

.0674 

.0808 

95.3 

07 

15 

23 

.4803 

-.1454 

.0419 

.0743 

91.7 

07 

16 

23 

.4933 

-.1553 

.0249 

.0714 

93.2 

07 

17 

23 

.5109 

-.1958 

-.0125 

.0829 

92.9 

07 

18 

23 

.5613 

-.2113 

-.0141 

.0909 

87.5 

07 

19 

23 

.5217 

-.2225 

-.0735 

.1073 

86.3 

09 

12 

23 

.5093 

-.2167 

.0188 

.0942 

91.3 

09 

13 

23 

.5061 

-.2078 

.0003 

.0832 

92.8 

09 

14 

23 

.4996 

-.1992 

-.0191 

.0873 

94.1 

09 

15 

23 

.4932 

-.2044 

-.0406 

.0980 

92.9 

09 

16 

23 

.5028 

-.2093 

-.4170 

.2505 

91.5 

07 

14 

16 

.5736 

-.0912 

.1236 

.0964 

93.6 

07 

14 

17 

.5907 

-.0966 

.0945 

.0800 

94.1 

07 

14 

18 

.5384 

-.1119 

.0787 

.0772 

93.7 

07 

14 

19 

.5523 

-.1169 

.0748 

.0762 

93.4 

07 

14 

20 

.4789 

-.1323 

.0622 

.0769 

94.4 

07 

14 

21 

.5186 

-.1273 

.0696 

.0776 

93.9 

07 

14 

22 

.4507 

-.1739 

.1090 

.1117 

94.8 

07 

14 

23 

.4726 

-.1387 

.0674 

.0808 

95.3 

07 

14 

24 

.5417 

-.1371 

.0726 

.0821 

95.6 

07 

14 

25 

.4528 

-.1491 

.0673 

.0847 

95.8 

07 

14 

26 

.4945 

-.1069 

.0258 

.0511 

96.0 

07 

14 

29 

.5009 

-.1409 

.0564 

.0758 

96.1 
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Table 3.1A (Continued) 
Sensitivity Testing Summary: M 


GOES Times Output 



(LST) 


Value 

1 

2 

3 

M-.5 

08 

11 

23 

.5596 

08 

12 

23 

.5090 

08 

13 

23 

.4951 

08 

14 

23 

.4972 

08 

15 

23 

.5021 

08 

16 

23 

.4802 

08 

17 

23 

.5190 

08 

18 

23 

.5230 


Worst 

Storage 

9M 3M Error 

3DTD 3DTN M»0.5 


-.2129 .0733 .1106 
-.1540 .0533 .0793 
-.1521 .0295 .0715 
-.1383 .0135 .0695 
-.1638 -.0082 .0684 
-.1875 -.0415 .0911 
-.2090 -.0435 .1003 
-.1708 -.0798 .0964 


2 

r 

Percent 


95.3 

94.9 

89.7 
90.2 

92.9 
93.0 

91.7 

90.9 
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Table 3. IB 


Sensitivity Testing Summary: P 

Worst 

GOES Times Output Storage 


1 

(LST) 

2 

3 

Value 

M-.5 

3M 

3DTD 

3M 

3DTN 

Error 

M-0.5 

t. 

r 

Percent 

06 

14 

23 

.06529 

.00589 

-.01140 

.00700 

99.3 

07 

14 

23 

.06709 

.00469 

-.00989 

.00584 

98.4 

08 

14 

23 

.06556 

.00472 

-.00824 

.00523 

98.1 

09 

14 

23 

.06549 

.00807 

-.00722 

.00723 

99.5 

10 

14 

23 

.06817 

.01034 

-.00618 

.00666 

99.7 

11 

14 

23 

.06699 

.01326 

-.00614 

.00781 

99.5 

07 

09 

23 

.06159 

-.00595 

-.00453 

.00389 

99.6 

07 

10 

23 

.06131 

-.00162 

-.00590 

.00286 

99.5 

07 

11 

23 

.06388 

.00121 

-.00673 

.00326 

99.5 

07 

12 

23 

.06252 

.00285 

-.00813 

.00445 

99.4 

07 

13 

23 

.06422 

.00398 

-.00915 

.00528 

99.1 

07 

14 

23 

.06709 

.00469 

-.00989 

.00584 

98.4 

07 

15 

23 

.06978 

.00579 

-.00996 

.00637 

97.1 

07 

16 

23 

.06946 

.00687 

-.01120 

.00731 

96.5 

07 

17 

23 

.06771 

.01140 

-.01314 

.01015 

98.4 

07 

18 

23 

.06861 

.01323 

-.01747 

.01274 

99.0 

07 

19 

23 

.06938 

.00833 

-.03391 

.01715 

98.5 

09 

12 

23 

.06188 

.00397 

-.00741 

.00460 

99.4 

09 

13 

23 

.06160 

.00648 

-.00765 

.00575 

99.4 

09 

14 

23 

.06549 

.00807 

-.00722 

.00623 

99.5 

09 

15 

23 

.06636 

.00941 

-.00750 

.00687 

99.6 

09 

16 

23 

.06574 

.00970 

-.00920 

.00767 

99.3 

07 

14 

16 

.05235 

-.00054 

.02752 

.01055 

95.7 

07 

14 

17 

.05787 

.00184 

-.01872 

.00863 

96.5 

07 

14 

18 

.05510 

.00290 

-.01473 

.00729 

97.2 

07 

14 

19 

.06200 

.00303 

-.01133 

.00586 

97.8 

07 

14 

20 

.06125 

.00311 

-.00966 

.00581 

98.1 

07 

14 

21 

.06503 

.00340 

-.00983 

.00535 

98.0 

07 

14 

22 

.06685 

.00423 

-.00992 

.00569 

98.3 

07 

14 

23 

.06709 

.00469 

-.00989 

.00584 

98.4 

07 

14 

24 

.06770 

.00499 

-.00958 

.00583 

98.4 

07 

14 

25 

.06805 

.00523 

-.00930 

.00581 

98.5 

07 

14 

26 

.06835 

.00539 

-.00904 

.00576 

98.6 

07 

14 

29 

.06791 

.00574 

-.00869 

.00575 

98.8 
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less-acurate GOES Imagery is used. Since the preferred testing 
strategy cannot be used, another strategy must be devised. 

Previous work (Carlson and Boland, 1978, Dodd, 1979) indicates 
that the numerical boundary layer model used in CD and the GOES 
method provides acceptable surface temperature values for a vide 
range of soil parameters. It was decided to use the numerical 
model to simulate temperature data for GOES method testing. Since 
M and P are input parameters of the model, GOES method testing was 
restricted to diagnosing M and P. A master run of the numerical 
model was executed to provide surface temperatures at 48 GOES times 
for 63 ordered pairs of (M, P) values. Meteorological and geographi- 
cal initializations were obtained from a St. Louis case studied by 
Carlson et al. (1981) and Kocin (1979). The output from the master 
run was used to represent GOES imagery of a fictitious target area 
exhibiting considerable variability in M and P. A test run was then 
executed, using the 16 (M, P) pairs to be used in the subsequent GOES 
method case studies (Section 4) . The M and P values in the test 
run differed from those used in the master run. GOES time triplets 
were selected, and DTD and DTN were computed from the output of the 
test run. Bicubic regression equations were produced for M and P 
from the output of the test run. These regression equations were 
applied to appropriate DTD and DTN values obtained from the master 
run output. The resulting M and P values were then compared to the 
M and P inputs used to generate the master run output. 
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3.6.2 Sensitivity Studies Results 

The results of the sensitivity and accuracy studies are presented 
in Table 3.1. 

Actual values of M and P from the test-data regression equations 
and master-data temperatures are presented in Column 4 of Table 3.1. 
Relatively little importance is attached to these figures. The 
temperatures passed to the regression-evaluation routine were not 
observed temperatures; they were model-output simulated temperatures 
from the master run. Actual data could produce different results, 
due to differing amounts and types of error in model output and 
GOES imagery. If the results tabulated in Column 4 were perfect, 
all one could conclude is that the regression routine represents 
the model perfectly at (M, P) • (0.5, 0.065). Although the ability 
of the GOES method to represent itself is highly desirable, that 
ability is net essential to the work at hand. The quantities 
tabulated in Column 5 and 6 are more important to the current study 
than the values tabulated in Column 4. 

Columns 5 and 6 contain the partial derivatives of the regres- 
sion equation with respect to the predictands DTD and DTN respectively. 
These partials were evaluated with simulated temperatures from the 
(M, P) = (0.5, 0.065) master model run. They represent the sensi- 
tivity of the GOES method to small perturbations (or errors) in 
model output and/or data. In a typical model run, when (M, P) varies 
from (0.05, 0.0125) to (0.95, 0.125), DTD and DTN vary by about 20 K. 

A value of Q DTp of 0.05/K or of 0.005/K is sufficient to account 
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for the real variation of H and P as functions of DTD and DTN. The 
tabulated values are generally larger, Indicating possible numerical 
difficulties. 

Column 7, labeled ’worst storage error,' contains 

°- 4 IjotoI + laEBil or °' 4 hiSffil + ■aslffl • Thl ‘ < >“ ntlty 18 

,iie maximum error incurred by storing GOES temperatures in single 
bytes of computer stor ge. The GOES imagery used in this thesis was 
obtained from th ‘ McIDAS facility at the University of Wisconsin 
in Madison, Wl. Each pixel of the satellite image was allocated 
one byte of magnetic tape storage. The raw intensity levels are 
converted into Kelvin temperatures by a preprocessing routine. The 
temperatures generally fall into a 100 K range: from 243 to 343 K. 

If the mapping from raw intensities to temperatures is nearly 
linear and nearly the entire domain of intensities is mapped into 
nearly the entire 243-343 K range, the temperatures can be no more 
accurate than + 0.2 K. The subtraction step to form DTD and DTN 
doubles this error range. If the numerical model and the satellite 
data were completely error-free, the GOES method could still have 
errors as large as the values in Column 7, due to the method of 
storing data on the McIDAS input tapes. 

2 

Column 8 displays the multiple regression coefficient, r ~ , 
a statistic generated by the MINITAB statistical package (Ryan, 

Joiner, and Ryan, 1976). It is a measure of the quality of the 

2 

regression equation, r is the square of the correlation coefficient 
between the M (o 1 a ) values in the input to MINITAB and the M (P) 
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values obtained by substituting the DTD and DTN values in the MINITAB 
input into the egression equation for M (P) produced by MINITAB. 

r"~ is the fraction of variance explained by the regression equation. 

2 

A value of r* close to 100 percent indicates that the regression 

equation can be used to invert the model without appreciable error. 

The results of the GOES times screening study were evaluated 

subjectively. Three objective criteria were applied, but these 

criteria were selected subjectively; 

2 

•avoid time triplets wehre r for M is less than 90 percent 

•avoid worst-storage-errors for M exceeding 0.10 

•avoid worst-storage errors for P exceeding 0.010 

2 

GOES time triplets were rejected when r“ began to decrease signifi- 
cantly as one of the times was varied away from its optimum. 

Triplets were rejected when sensitivities (Column 5, 6, and 7) 
began to increase rapidly as one GOES time was changed . Also con- 
sidered, but not shown in Table 3.1 was the performance of the 
method at extreme values of (M, P) . Pronounced difficulties were 
occasionally encountered for (M 0.8; ! > 0.10), where the inversion 
of the model is numerically least stable. 

The results of the sensitivity study are summarized below: 

1) The preferred midday GOES time is after local noon, 
preferably between 1300 and 1500 but no later than 1600 
local standard time (1ST) . The time of surface tempera- 
ture maximum is optimum. 

2) Acceptable morning GOES times are after sunrise, but before 
,1000 1ST. The optimum is about 1 1/2 h after sunrise. 
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Times after 0900 will provide marginally acceptible results, 
if afternoon and evening times are optimal. 

3) The night GOES time should be at least 2 h after sunset. 
Performance of the method improve s slowly as the night GOES 
time is made later than 2300 LST. Model execution time 
(hence, cost) also increases as the night time is made later. 
The marginal Improvement in performance probably does not 
justify the expense of selecting night images after local 
midnight in an operational setting. 

4) If the midday image time deviates more than 30 minutes 
from optimum, flexibility in the choice of morning image 
time is reduced. 

These results were obtained from model simulations performed 
on the St. Louis case used in Carlson et al. (1981) and Kocin (1979). 
They may not have general applicability. Rose (1982, personal 
communication) reports poor results with the Kansas case study 
(Section 4), using GOES times of 0830, 1400, and 0100 1ST. The 
sensitivity studies of the St. Louis simulation indicate that (0830, 
1400, 0100) should be suboptima' but acceptible. This failure 
may be related to the choice of M and P values used in the model. 

Rose is using a smaller P range and a more sophisticated algorithm 
for generating (M, P) pairs for model initialization. Additional 
experimentation with the GOES method should lead to refinements in 


the criteria cited above. 
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Another restriction on the choice of GOES times was discovered 
during the processing of the Kansas case study (Section 4) : 

5) If possible, the model should be initialized to start after 
effective sunrise, about 1 h after astronomical sunrise. 

The morning GOES time should be at least 30 minutes after 
model- initialize ion time. If this is impossible (due to 
limited imagery or adverse cloud conditions), then: 

a. Initialize the model to start at the selected GOES 
time 

b. Use the model initial conditions (the initial surface 
temperature, etc.) as the model output for the first 
GOES time 

Effective sunrise of radiation sunrise is defined as the time at 

which the surface first receives an excess of input energy (incident 

radiation plus G ) over emitted radiation. At effective sunrise, 
o 

H q becomes positive. If the model is initialzed to start before 
effective sunrise, the model enters the nocturnal mode and produces 
a very shallow, spurious nocturnal inversion. Modeled surface 
temperatures are unreliable until this inversion is destroyed by 
heating. Use of model initial conditions as model output is not 
desirable, but it is better than the use of purely erroneous data 
resulting from the nocturnal-mode computations encountered if the 
model is started before effective sunrise. 

For a few GOES time triplets in the sensitivity study, MINITAB 


produced the diagnostic message: VARIABLE . . . OMITTED FROM 
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REGRESSION EQUATION. This message implies that the model output 

for that GOES time triplet is overfit by a bicubic. One of the 

2 

cubic terms is discarded by MINITAB. When this occurs, r decreases 
and the sensitivity measures (Column 5, 6, 7) increase. There is no 
a priori way of determining which GOES time triplets will suffer 
from this overfitting. It has not been observed to happen in an 
optimal GOES time triplet; it occurs infrequently in suboptimal but 
acceptable triplets. The problem is almost certainly situation- 
dependent: different geographical and/or meteorological 

initializations imply different GOES time triplets having this 
property. MINITAB is most likely not in error. By coincidence, 
the model output has essentially no cubic component in one of the 
predictors, DTD or DTN. To eliminate this ca*,se of impaired per- 
formance, possible GOES triplets could be prescreened by the routine 
used in the sensitivity study: 

1) Initialize the model with the appropriate geographical 
and meteorological parameters, including the M and P 
values to be used in the actual operational run. 

2) Prepare a test dataset spanning all GOES times which are 
available for use. 

3) Select a GOES time triplet from the available GOES times, 
subject to the constraints mentioned above in this 
section. Select the most nearly -optimal first. 
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4) Compute regression equations for any predictand, for the 
GOES triplet in question. 

5) If the 'VARIABLE . . . OMITTED . . message appears, 
return to Step 3. Otherwise, accept the GOES time triplet. 
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4.0 RESULTS OF TEST CASES 
4.1 INTRODUCTORY REMARKS 

This thesis is an attempt to answer two major questions: 

1) Can M, P, H q , and/or E be diagnosed on the regional scale 
via satellite imagery? 

2) Can GOES imagery be incorporated into the CD method? 
Preliminary work (Section 3) indicates that both questions can 
be answered in the affirmative. It is necessary, however, to 
demonstrate the GOES method on actual regional data before making 
any meaningful claims. 

Two test cases were considered: 

1) Indiana/Illinois/Kentucky , 22 August 1978 

2) Eastern Kansas, 27 July 1978 

These cases were selected because both GOES and HCMM data were 
available. The Kansas case was also selected as a pilot for a 
future study of time-variability of M and P anomalies (Rose, 1983). 
Since none of the diagnosed variables could be measured on the 
regional scale, verification of the results of the test cases is 
difficult. The output of the two competing methods, CD and the 
GOES method, will be presented, along with crop-moisture-index and 
the preceding three week’s rainfall. To the extent that both methods 
agree and are consistent with the available ground measurements, 
the GOES method can be judged a useful method. 
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The Indiana/Illinois/Kentucky case (henceforth called the 
Indiana case) is discussed in Section 4.2; the Kansas case, in 
Section 4.3. Figures 4.1 through 4. 25, illustrating the Indiana 
and Kansas cases, are collected in Section 4.4 for convenience. 
General conclusions are presented in Section 5.1. 

4.2 THE INDIANA CASE: 22 AUGUST 1978 

So many things I could have done/But clouds got in my way 

Judy Collins 

Showers and thunderstorms left many bands and patches of 
cloud over the Indiana study area. At four of the five image times 
used in the study, there were clouds obscuring significant portions 
of the study area. The only clear image was the HCMM night image, 
0748 Z, 21 August 1978, 48 hours before the modeled image time. 

Under these circumstances, the extra image used in the GOES method 
is a disadvantage. Whenever clouds appear on any image used, output 
values are meaningless. The third GOES image introduces a third 
set of clouds into the diagnostic method, obscuring more of the 
region. Only about 25 percent of the GOES/Indiana domain was not 
obscured by heavy clouds; thin clouds may affect about 10 percent 
more of the domain. About 50 percent of the HCMM/Indiana domain is 
usable. 

Data results for the Indiana case are presented in Figures 4.1 
through 4.11. There is evidence that both methods successfully 
diagnosed large-scale moisture-availability patterns, despite the 
interference by clouds. The three-week rainfall charts (Figures 4.2 
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and 4.7) show a marked antecedent precipitation minimum near the 
intersection of the Wabash and Ohio Rivers. Three-week precipita- 
tion and crop-moisture index* both indicate greater moisture to the 
north and east; the variation in crop moisture index values is 
less dramatic than that of the antecedent precipitation. Figures 
4.4 and 4.9 show M as diagnosed from HCMM and GOES data, respec- 
tively. Both have dry values in the southwest quadrant. Neither 
is quantitatively accurate: M of 0.25 would indicate wilting 

vegetation, while M of 1.00 represents a wet surface or open water. 
The pattern of values correlates with the antecedent precipitation 
adequately, considering the masking due to clouds. 

The diagnosis of P, H q , and E in the Indiana case will not be 
discussed. The reader may discern for himself the degree of cor- 
relation between HCMM and GOES fields, and their relationship to 
three-week precipitation, crop moisture, and/or land-use patterns. 

P, H , and E will be discussed in the Kansas case (Section 4.3), 
o 

where interference by cloudiness is markedly less evident. 

The HCMM E field (Figure 4.6) contains a curious and counter- 
intuitive feature which will be discussed. A band of minimum E 
almost exactly parallels the Ohio River. One might expect high 
evaporation rates over bodies of water, all other conditions 
being equal. A definitive explanation for this minimum in E is 
not available. The river was probably cooler than the surrounding 

"^Operational product of NOAA/USDA Joint Agricultural Research 
Facility; broadcast as NAFAX map N118 and DIFAX map D262. 
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land during the daylight hours, although the satellite imagery 
failed to resolve a temperature anomaly. Evaporation from cool 
water could be lower than evaporation from warm land. The night 
HCMM image indicates a temperature maximum in the Ohio River Valley. 
This temperature maximum apparently was transformed into the E 
minimum by the regression routine. If total evaporation was actually 
lower over the river than over the land due to cool daytime water 
temperatures, the model has obtained a correct result for the wrong 
reason. If the evaporation minimum actually did not exist, the 
model has generated a spurious feature, 

4.3 THE KANSAS CASE: 27 JULY 1978 

The Kansas data are much better than the Indiana data. Clouds 
obscured less than 25 percent of the image area. As a result, 
all output fields could be obtained and compared. Data and results 
for the Kansas case are presented in Figures 4.12 to 4.25. 

Meteorological conditions over the eastern Kansas study area 
were favorable. At 2100 Z, 27 July 1978, winds were generally 
5-10 knots from the north to northeast. The pressure gradient was 
weak. By 0900 Z, 28 July 1978, winds had shifted to southerly, 
but had decreased to 0-5 knots over the area of interest. The 
temperature and wind patterns shown in Figure 4.13 are consistent 
with an approximate nonadvective assumption. Cloud conditions are 
not plotted in Figure 4.13 but satellite Imagery shows relatively 
few patches of clouds over the study area. The three— week antecedent 
precipitation maps (Figures 4.14 and 4.20) indicate a dry region in 
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the southwest quadrant of the study area. Crop moisture Index shows 
a similar dry area, displaced a bit eastward. A region of high 
crop moisture is found in the northwest quadrant of the region, 
consistent with the three-week precipitation but not entirely 
explained by it. 

Several differences between HCMM and GOES temperature fields 
are evident in Kansas case Imagery. The superior resolution of the 
HCMM satellite is clearly demonstrated. None of the GOES Images 
(Figure 4.21A-C) displays any temperature feature which can be 
identified with any of the lakes in the study area. Both HCMM 
images (Figure 4.15A-B) show one or more lakes. Note in particular 
Perry Reservoir, the lake in the northeast quadrant of both the 
HCMM and the GOES domains. It appears as a distinct temperature 
minimum on the 1934 Z HCMM image and as a distinct maximum on the 
0823 Z HCMM image. The high effective P of open water implies that 
lakes should be cooler than the surrounding land during the day 
and warmer than the land at night. For small bodies of water HCMM 
is able to observe this distinction; GOES is not. The 8 km resolution 
of GOES may explain its inability to see the lakes; Perry Reservoir, 
however, is larger than 8 km in dimension. 

The temperature patterns displayed on the 1934 Z HCMM image 
(Figure 4.15A) and the 2000 Z GOES image (Figure 4.21B) are similar, 
but the HCMM image is about 10 C (10 K) warmer than the GOES image. 

The HCMM image was collected on 2,8 July 1978; the GOES-, on 27 July 
1978. A 10 C change in surface temperature between the two image 
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times seems excessive, considering the fairly stable weather pat- 
tern over the region. Both images are considerably warmer than the 
surface air temperatures displayed in Figure 4.13A. This is to be 
expected: the actual surface of the earth, exposed to direct sun- 

light, should be considerably warmer than the air inside an 
instrument shelter 1.3 m above the heated surface. A difference 
of about 25 C, as implied by the HCMM imagery, is possible; 
intuitively, 25 C seems extreme. 

The land-use patterns (Figure 4.12) are not well correlated 
with the temperature fields. Most of the region consists of some 
mix of unirrigated cropland and rangeland. The principal crop in 
eastern Kansas is wheat. Wheat is a grass. The distinction between 
wheat and range grass is obvious to the farmer and the cartographer; 
the distinction may be much less obvious to the longwave-IR sensor 
of a weather satellite. Unfortunately, there was no significant 
area of irrigated cropland in the study area. Irrigated cropland 
would be expected to differ in surface temperature from unirrlgated 
land. 

The four output fields show good agreement between HCMM and 
GOES methods, good consistency among field types, and good agreement 
with three-week antecedent precipitation. The agreement with crop 
moisture index is adequate, but not as good as the relationship 
with antecedent rainfall. The output fields are not to be inter- 
preted as exact values, though. The range of M displayed in 
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Figures 4.16 and 4.22 is unrealistically large, for example. This 
qualitative accuracy and quantitative inadequacy also occurred 
in the Indiana case (4.1). 

Both M fields clearly show the dry region in the southwest 
of the HCMM domain and in the west of the GOES domain. Clouds 
unfortunately obscure the region of highest M gradient in the HCMM 
domain. The HCMM M field (Figure 4.16) resolves the Perry Reservoir 
asanM maximum; the lake on the northwest edge of the domain may 
be resolved. 

The P fields show low values of P in the dry region, as would 
be expected. The HCMM P field (Figure 4.17) resolves two lakes as 
P maxima; the GOES field (Figure 4.23) does not resolve the lakes. 

The H q fields (Figures 4.18 and 4.24) show high heafi flux over 
dry terrain. 

The E fields (Figures 4.19 and 4.25) show the expected result: 
low evaporation in the dry region. Neither HCMM nor GOES E field 
resolves either lake. The E fields correlate better with previous 
three-week rainfall than with crop moisture index. 

The reader is encouraged to study the output fields, the input 
fields, and the meteorological and physiographic data and reach 
his own conclusions about the relative merits of the GOES method 
and the CD method for diagnosing M, P, H q and/or E. 
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4.4 FIGURES DEPICTING CASE STUDY RESULTS 

Figures 4.1 through 4.25 uepict initial conditions and results 
of the two case studies described in Sections 4.2 and 4.3. The 
figures are presented in the following order: 

Indiana case: Figures 4.1 through 4.11 

4.1: Base map and land-use data 

4.2 through 4.6: Indiana-case HCMM: 

4.2: Three-week rainfall and crop moisture index 

4.3: Satellite images (temperature fields) 

4.4 through 4.6: Output fields (M, P, E) 

4.7 through 4.11: Indiana-case GOES: 

4.7: Three-week rainfall and crop moisture index 

4.8: Satellite Images 

4.9 through 4.11: Output fields (M, P, E) 

Kansas case: Figures 4.12 through 4.25 

4.12: Base map and land-use data 

4.13: Surface weather analyses 

4.14 through 4.19: Kansas-case HCMM: 

4.14: Three-week rainfall and crop moisture index 

4.15: Satellite images 

4.16 through 4.19: Output fields (M, P, H^, E) 

4.20 through 4.25: Kansas-case GOES: 

4.20: Three-week rainfall and crop moisture index 

4.21: Satellite images 

4.22 through 4.25: Output fields (M, P, H q , E) 



On all satellite imagery and output fields, 
by cloud cover are indicated by scalloped lines; 


areas affected 
contours are 


omitted in cloudy areas. 
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Figure 4.1A: Geographical location of GOES and HCMM domains for 

the Indiana test case. 


Heavy solid lines: GOES domain 

Heavy dashed lines: HCMM domain 

Light solid lines: Rivers and lake shores 

Light dashed lines: State boundaries 
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Figure 4. IB: Land use map for th< 

A: Cropland and woodland 1 

B: Cropland 1 
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Figure 4.2: Three-week antecedent rainfall and crop moisture index 

for the Indiana case HCMM domain. 


Solid contours: Preceding three-week precipitation (inches) 

Dashed contours: Crop moisture index 
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Indiana-case HCMM image, 1858Z, 22 August, 1978. 
Contours are satellite-derived surface temperatures, 
in °C, corrected for atmospheric water vapor. 
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Figure 4.3B: Indiana-case HCMM image, 0748Z, 21 August, 1978. 

Contours as in Figure 4.3A. 
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Figur e 4.4: Indian.' i-case noisture availability diagnosed iron HCMM 

S data. Contour interval: 0.25. 


52 


ORIGINAL. PAGE 19 

OF POOR QUALITY 



r — r "T — i 

0 100 KM 

Figure 4.5: Indiana-case thermal inertia (P) diagnosed from HCMM 

data. Contour interval: 0.02 cal cm K sec - !' 2. 
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Figure 4.6: Indiana-case total evaporation diagnosed from HCMM 

data. 
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Three-week antecedent precipitation and crop 
moisture index for the Indiana case GOES domain. 
Contours as in Fig. 4.2. 


0 

Figure 4.7: 
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Figure 4.8A: Indiana-case GOES image, 1230Z, 22 August 1978. 


Contours as in Fig. 4.3A. 
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Figure 4.8B: Indiana-case GOES image, 1900Z, 22 August 1978. 

Contours as in Fig. 4.3A. 
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Figure 4.8C: Indiana-case GOES image, 0600Z, 22 August 1978. 

Contours as in Fig. 4.3A, f 
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Figure 4.9: Indiana-case moisture availability <M) diagnosed from 

GOES data. Contours as in Fig. 4.4. 
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Figure 4.10: Indiana-case thermal inertia (P) diagnosed from GOES 


data. Contours as in Fig. 4.5. 
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Figure 4.11: Indiana-case total evaporation (E) diagnosed from 

GOES data. Contours as in Fig. 4.6. 
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Reference: Kansas Land Use Patterns , Summer 1973. Donald L. Williams 

and Bonnie L. Barker. Space Technology Laboratories , 
University of Kansas in Cooperation with Planning Division 
Kansas Department of Econ. Development. 
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Figure 4.13A: Surface weather analysis. Contours: pressure in mb. and 

tenths; thousands and hundreds digits omitted. 
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Figure 4.14: Three-week antecedent rainfall and crop moisture 

index for the Kansas case HCMM domain. Contours 
as in Fig. 4.2. 
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Figure 4.15A: Kansas-case HCMM image, 1934Z, 28 July 1978. 

Contours as in Fig. 4.3A. 
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Figure 4.15B: Kansas-case HCMM image, 0823Z, 27 July 1978. 

Contours as in Fig. 4.3A. 
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Figure 4.16: Kansas-case moisture availability (M) diagnosed from 

HCMM data. Contours as in Fig. 4.4. 
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Figure 4.18: Kansas-case heat flux (H ) diagnosed from HCMM data. 

Contours in W m“ . 

1: 0.15 2: 0.20 3: 0.25 4: 0.30 

5: 0.35 6: 0.40 7: 0.45 
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Figure 4.19: Kansas-case total evaporation (E) diagnosed from 

HCMM data. A special model run using 32 (M, P) 
pairs was used to produce this diagnosis. 
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Figure 4.20: Three-week antecedent rainfall and crop moisture index 

for the Kansas case GOGS domain. Contours as in 
Fig. 4.2. 
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Figure 4.21A: Kansas-case GOES image, 1100Z, 27 July 1978. 

Contours as in Fig. 4.3A. 
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Figure 4.21B: Kansas-case GOES image, 2Q00Z, 27 July 1978. 

Contours as in Fig. 4.3A. 
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Figure 4.21C: Kansas-case GOES image, 0900Z, 28 July 1978. 

Contours as in Fig. 4.3A. 
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Figure 4.22: Kansas-case moisture availability (M) diagnosed d 


from GOES data. Contours a? in Fig. 4.4. 
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0 100 KM 

Figure 4.23: Kansas-case thermal inertia (P) diagnosed from GOES 


data. Contours as in Fig. 4.5. 
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0 100 KM 


Figure 4.24: Kansas-case heat flux (H q ) diagnosed from GOES data. 

Contours in W m"'* ; 

0: 0.05 1: 0.15 2: 0.25 3: 0.35 4: 0.45 

Note that contour 1 of Fig. 4.18 corresponds to contour 1 of this 
figure; contour 7 of Fig. 4.18 corresponds to contour 4 of this 
figure. 
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0 100 KM 

Figure 4.25: Kansas-case total evaporation (E) diagnosed from 


GOES data. Contour interval: 0.10 cm. 
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5.0 CONCLUSIONS AND SUGGESTIONS FOR FUTURE RESEARCH 
5.1 CONCLUSIONS 

The following conclusions can be drawn: 

1) Both the GOES method and CD can be used for regional-scale 
studies 

2) Three GOES images are required for the GOES method. The 
optimal image times are: 

a) about 90 min after sunrise 

b) the time of maximum surface temperature 

c) around midnight 

There is some flexibility in the choice of image times , 
particularly in the choice of the third image time. 

3) The superior resolution of HCMM is important, even on the 
regional scale. Significant features can be resovled via 
CD which are invisible to the GOES method. 

4) The superior resolution of HCMM enables CD with only two 
images to equal or exceed the performance of the GOES 
method with three images. 

5) Relative diagnoses of M, P, and H q can be obtained from 
either method. The ability to diagnose E has not been 
satisfactorily demonstrated. Although the E fields 
produced by both methods appear qualitatively reasonable 
in most areas and are probably adequate, there are some 
inconsistencies that need to be explained. 
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6) Neither method may produce quantitatively correct output 
fields. The diagnosed M fields exaggerate naturally* 
occurring contrasts. The same type of error probably 
occurs in the E and H fields. Analysis of error in the 
P fields is complicated by the difficulty of defining P 
in a vegetative canopy. Negative and near-zero values 
of P are certainly incorrect. 

7) The present work must be considered a preliminary study. 

The development of the GOES method is far from complete. 

5.2 LIMITATIONS OF THE METHOD: DISCUSSION AND SUGGESTIONS FOR 

FUTURE RESEARCH 

Considerable research is needed in defining and measuring 
surface parameters in vegetative canopies. Diagnoses of these 
parameters must remain tentative until more is known about them. 

The difficulties encountered in defining and treasuring surface 
temperature have been discussed (Section 3.4). Several other 
questions involving vegetative canopies must be addressed: 

1) How is P defined in a vegetative canopy where the radiating 
surface is a significant distance above the solid ground 
surface? How is the temperature of the radiating surface 
influenced by heat transport to and from the region below? 

2) To what extent is the temperature of the effective radiating 
surface affected by ventilation, when the radiating surface 
is at some distance above the ground in a canopy? 

3) Does M vary significantly during the day as plant stomata 
open and close in response to a variety of stresses? 
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4) Does E vary significantly as plant stomata open and close? 

Is any significant part of such variation independent of 
the variation of M? 

These and similar questions need to be answered through further 
observations in vegetated areas. Significant input into the 
research effort must be provided by plant physiologists. The 
research will be neither easy nor inexpensive. 

To date, neither CD nor the GOES method has been attempted on 
an autumn or winter case. The numerical considerations cited in Sec- 
tion 3.4 suggest that some difficulty may be encountered in diagnosing 
the surface parameters from autumn or winter imagery, even if the 
question of the thermal and radiative properties of snow is ignored. 
The daily range of surface temperature may be reduced to the point 
that the errors inherent in the model overwhelm the information in 
the satellite images. The numerical model has been demonstrated 
successfully on simulated autumn initial conditions. The question 
of whether regression equations obtained from such initializations 
are adequately stable to permit the mapping of satellite images 
into M, P, H q , and/or E fields has not yet been addressed. 

The superior resolution of HCMM proved to be of value even 
in the regional-scale studies (Section 4) . The ability to resolve 
a feature like the Ohio River is certainly worthwhile. HCMM is no 
longer operating, however, and there is no likelihood that a 
similar mission will be launched in the near future. The CD 
algorithm could be tested on TIROS-N imagery or on other thermal-IR 
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imagery from polar-orbiting platforms. Although TIROS-N has poorer 
resolution than HCMM, its resolution (1 km) is distinctly superior 
to the resolution of GOES. The question of whether the superior 
resolution adequately compensates for the paucity of overpass 
times is one that can only be answered by experiment. 

Hybrid methods could also be attempted to take advantage of 
the availability of GOES imagery while retaining the superior 
resolution of polar-orbiter Imagery. On a day when a polar-orbiter 
image pair is available, for example, one GOES image could be used — 
probably the early-morning image — to provide the three images 
needed for the GOES method. On a day with a single polar-orbiter 
overpass, two GOES images could be combined with the single polar- 
orbiter image. It is not known whether differences in image quality 
or calibration are large enough to prevent the mixing of image 
sources in the GOES method, nor is it known whether the difference 
in resolution between GOES and polar-orbiter imagery will make the 
two types incompatible. 

The image differencing done in the GOES method has a theoretical 
deficiency which was not realized until after the completion of the 
case studies presented in Section 4. Image subtraction is a valid 
technique only under the conditions cited in Section 3.5. The 
assumption of similar soundings at all points In the domain is 
crucial. If clear skies and light winds prevail during the night, 
this assumption is likely to fail. Cold-air drainage will produce 
local pockets of anomalously cold air at the surface; just above 
these shallow cold layers, there will be no horizontal contract. 
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These shallow cold pockets are definitely large enough to be 
resolved by HCMM; some are resolvable by GOES.* This error in 
the subtraction procedure is maximized by selecting a morning image 
time near sunrise; the error decreases as the morning image time 
is made later. Unfortunately , a later morning image time increases 
the numerical difficulties cited in Section 3. A. The trade-off 
between inaccuracies and numerical difficulties could be investi- 
gated, but only if some form of ground truth is available. 

Preliminary evidence indicates that the numerical difficulties are 
the more serious. There is no evidence of degradation of performance 
attributable to cold pockets in either of the case studies in 
Section A. 

The subtraction step could be bypassed, and ull three tempera- 
tures could be used as predictors in a modified GOES procedure. If 
this were done- the number of coefficients in the regression 
equations would increase, increasing the risk of overfitting the 
data (Section 3. A and Appendix X). A.i&o, the utt&ssflying problem 
would not be addressed. The riddel would still be initialized with 
only one sounding. The method would implicitly assume that the 


An extreme example of this phenomenon is presented by Schlegel and 
Butch (1980) . They found that the nocturnal temperature drop is 
greatly exaggerated in the Barrens, a notorious cold-air drainage 
location near State College, PA. Minimum temperatures in the 
Barrens are much lower than minima at Penn State University on 
clear, calm nights. Daytime maxima in the two locations are com- 
parable. Much of this anomaly is explained by cold-air drainage 
into the bowl- shaped valley. The thermal properties of the sandy 
soil and stunted vegetation cannot alone account for an effect of 
the magnitude of the cold anomaly of the Barrens. 
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sounding was identical at all points in the region: a more stringent 

assumption than the assumption of similar— shaped soundings. However , 
the quality of the output might bu Improved. A test of a three- 
temperature GOES method is desirable, but should not be given high 
priority. Unless a new idea is forthcoming, the error in the sub- 
traction algorithm will probably have to be accepted. 

The one-dimensional boundary-layer model used, in CD and the 
GOES method has proved to be adequate in diagnosing M, P, H q , and E 
from HCMM and GOES imagery in several different studies. Work is 
presently in progress to add several improvements to the model: 

1) constant temperature advection, 

2) improved daytime boundary-layer winds, 

3) better subsurface heat storage and flux parameterization. 
These improvements should be documented by Rose in 1983. 

Both CD and the GOES method use a simple polynomial regression 
equation in two variables to invert the model. Results have been 
qualitatively satisfying but quantitatively dubious. This 
inadequacy may be related to the physical response of surface 
temperature to M and P — a response which is apparently modeled 
quite well. 

Diurnal surface temperature rise and fall are extremely sensi- 
tive to small changes of M and P when M and P both have relatively 
low values. For high values of M and P, diurnal temperature rise 
and fall are rather insensitive to changes in M and P. Note the 
data in Table 5.1, taken from a model simulation of the Kansas case 
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Table 5.1 

Model surface temperatures generated for the Kansas case 
(27 July 1978) for selected initial values of M and P 


SIMULATED SFC TEMPERATURES (K) 
AT SIMULATED GOES TIMES LST 


M 

P 

0800 

1400 

2300 

DTD 

DTN 

0.05 

0.0125 

303.21 

332.37 

292.89 

29.16 

39.50 

0.35 

0.0125 

299.61 

317.39 

289.45 

17.78 

27.94 

0.05 

0.0500 

296.62 

322,52 

300.18 

25.90 

22.34 

0.35 

0.0500 

295.75 

313.87 

296.57 

18.12 

17.3 

0.65 

0.0875 

294.10 

308.74 

296.91 

14.64 

11.83 

0.95 

0.0875 

293.75 

307.59 

296.66 

13.84 

10.93 

0.65 

0.1250 

293.89 

307.07 

297.50 

13.18 

9.57 

0.95 

0.1250 

293.59 

306.16 

297.28 

12.60 

8.91 


(4.2). When the model is inverted, the sensitivities are also 

inverted. Let AT be an appropriate measure of diurnal temperature 

9M 3P 

rise or fall. Then, and are very small for large AT, and 
rather large for small values of AT. It is difficult to obtain a 
single polynomial in temperature or temperature-change which matches 
the behavior of M or P at both ends of the temperature or temperature- 
change range. A piecewise-polynomial interpolating function may fit 
the data more successfully than a global least-squares approximating 
polynomial. The piecewise interpolant need not match the model- 
output globally. For example, a cubic spline interpolant can fit 
data very well, providing only a single predictor is used. If M 
were to be diagnosed from morning temperature change only, with 




P-dependence suppressed in Che model, a cubic spline could be used 
in place of: a regression equation in Che single predlctund DTD ? The 
extension of che concept of cubic splines to two dimensions is 
beyond Che competence of the present author. However, it might be 
possible to fit an augmented model-output data set by pieeewise- 
1 inear functions (see Appendix II). More than 16 (M, P) pairs would 
probabLy have to be used in the model simulation. The problem of 
replacing the regression routine by piecewise-poiynomial interpolates 
is worth pursuing. See Appendix II, Section 3 for a possible 
plecewisa-linear interpolation method. 

The chief impediment to the operational use of Cl) or the GOES 
method is the acquisition and alignment of satellite subimages. To 
obtain one aligned subimage of a target area, a technician must 
spend one half to one full working day searching for recognizable 
geographic features on the satellite image. The task is excep- 
tionally difficult on night imagery , where the contrast between 
land and ocean is minimal. Until the procedure for obtaining 
aligned subimages is expedited, real-time applications will be 
nearly impossible, and ongoing studies of several regions will be 
economically infeasible. High priority should be placed on 
developing a fully-automated procedure for selecting and aligning 
subimages for GOES and other infrared satellite, imagery. 

Although the GOES method is not yet ready for real-time 
operational use, work is currently being done on the time evolu- 
tion of the historical M anomaly in the Kansas case (Section 4.3) 
(Rose, 1983) . Other potential applications of the GOES method 
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Include the preparation of M, P, and possibly and E fields for 
use. by raesoBcale numerical weather prediction models. If the sub- 
image acquisition problem can be solved, the GOES method could 
provide surface parameters for operational NW1’ ; a pilot study should 
be done. 

5.3 FUTURE RESEARCH : A PLAN OF ACTION 

Remote sensing of soil and surface parameters is potentially 
very valuable. The GOES method has been shown to diagnose several 
soil and surface parameters with some skill, using weather-satellite 
imagery. The future usefulness of this technique depends on the 
results of further research. This author recommends pursuing the 
following course : 

As the top priority, solve the subimage acquis itiou and align- 
ment problem. Operational work cannot proceed as long as acquisition 
of usable data is extremely tedious and difficult, 

If funding and outside expertise can be found, the study of 
the mierometeorology of plant canopies should also be given high 
priority. It is hard to model a quantity which is poorly understood, 
using data of unknown quality and a method designed for conditions 
different from the ones which actually prevail. 

Applications of the GOES method and pilot studies of its uses 
should be given intermediate priority . 

Intermediate priority also should be given to improving the 
regression algorithm or replacing it with piecowise-polynomial 


functions , 
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Studies of TIROS-N and hybrid methods should receive high 
priority when the image acquisition problem is solved. The use 
of TIROS-N data in the CD method should be attempted as soon as 
data become available. Large-scale evaluation of TIROS-N and 
hybrid methods must wait for the solution to the image acquisition 
and alignment problem, for economic reasons. 

Low priority should be given to the incorporation of improve- 
ments into the model and the systematic evaluation of the image- 
differencing technique. Until the micrometeorology of plant 
canopies is understood, it will be difficult to evaluate the 
results of "improving" the model or to evaluate image-differencing. 
Preliminary experiments with "improvements" in the model indicate 
that uncertainties in the GOES method overwhelm any improvements 
resulting from added sophistication in the numerical oicdel. The 
time for major improvements in the model is after ground-truth 
becomes available, not before. 

There is little doubt that surface parameters will eventually 
be diagnosed via some form of remote sensing. Whether operational 
diagnoses can be achieved by the GOES method, using geostationary 
observation platforms, remains to be seen. Whether the GOES method 
or one of its successors is eventually adopted, the experience gained 
in the development of CD and the GOES method will be valuable. 
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APPENDIX I: NUMERICAL ANALYSIS 

Much of the work in this thesis involves approximations. The 
satellite data and model output are approximate; they contain error. 
This appendix discusses: 

1) theoretical implications of approximation 

2) the role of approximation in the GOES method 

3) possible effects of approximation 

A formal numerical analysis of the GOES method is not attempted. 
Basic principles of numerical analysis will be presented and 
applied to the GOES method. Awareness of these principles can 
suggest ways of reducing numerical errors in the GOES method. 

The important topic of floating-point arithmetic will not be 
discussed, since the numerical model used in the GOES method is 
relatively simple and stable. The interested reader is referred 
to Chapter 2 of Forsythe, Malcolm, and Moler (1977), Section 1.1 
considers the problem of approximating derivatives, and applies its 
results to the problem of approximating temperature derivatives 
with GOES data. Section 1.2 deals with regression: approximating 

one function by a simpler one, when only limited data of uncertain 
quality are available. 

1.1 NUMERICAL APPROXIMATION OF DERIVATIVES 

Assume that it is desirable to estimate the time derivative 
of temperature at a specific time: 0945 local time, for example. 
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The simplest and best approach would be to use the centered- 
difference approximation to the derivative: 


3T 

3T 


T 2 - T 1 
C 2 - h 


( 1 . 1 ) 


where and are temperatures measured at times and t 2 » 
respectively. Times t^ and are selected such that: 


2 (t l + t 2 ) = 0945 ’ 


( 1 . 2 ) 


This approximation is equivalent to approximating the graph of 
T(t) from t^ to by a straight line. 

Assume that the actual temperatures are as follows: 

Time Temperature 


0930 

1000 

The actual value of the derivative: 


300.0 K 

301.0 K 


(1.3) 


= 2.00 K/hr (1.4) 

Now, assume that satellite imagery or computer output is 
available, but that either contain random error of up to 0.3 K. 

This corresponds to a 0.1 percent error in absolute temperature, not 
unreasonable for a temperature computed from narrow-band IR reception 
by a satellite. The available data are: 

Time Temperature Range 

0930 299.7 K< •£ < 300.3 K (1.5) 

300.7 K <: T 2 £ 301.3 K 


1000 
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The derivative computed from this data could fall anywhere in the 
range: 0.8 K/hr to 3.2 K/hr. An estimate that crude would have 

little or no practical value. 

Consider the following data, containing the same uncertainty: 
Time Temperature Range (K) 

0730 295.7 < T, < 296.3 (1.6) 


204.7 < T 2 < 305.3 


1200 

The computed value of falls in the range: 1.87 K/hr to 2.31 K/hr. 

This estimate is useful. 

The above is an example of the problem of subtraction of 
nearly-equal approximate quantities. The relative error in the 
quantities can be greatly magnified. Finite-difference approxi- 
mation of derivatives is a numerically sensitive process because 
of this property of subtraction. 

An implicit assumption was made in the discussion of the second 
set of data: we assumed that the time derivative of temperature 

remained essentially constant for several hours. There is no 
a priori justification for that assumption; in fact, it is unlikely 

to be strictly true. Conte and De Boor (1972) demonstrate that the 

9T 

centered-dif ference approximation to has an error: 

2 


e- (-£> A 

6 9t 2 3 


) 


(1.7) 


t=t* 


where 


1 1 2 ~ t V = 2h 


( 1 . 8 ) 
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t* is some (unknown) time in the interval (t^, t^) , and T(t) 
must have three time derivatives throughout the interval. The 
error Increases as the square of the time interval between observa- 
tions. 

There is a trade-off between the accuracy of the centered- 
difference formulation and the numerical error Introduced by sub- 
traction of nearly-equal approximate quantities. In the situation 
encountered in this thesis, the errors from the subtraction are 
the more serious. 

1.2 REGRESSION, CURVE-ITITING, AND APPROXIMATE DAT/' 

Curve-fitting is a common problem in science. An experiment 
or simulation provides data, and the investigator wishes to express 
the data as a mathematical function. However: 

1) The data are not available for all possible experimental 
conditions. 

2) The data invariably contain some error or uncertainty. 

The investigator would like to obtain a mathematical result 

of general applicability which is consistent with the experimental 
or simuluated data. She is willing to accept a function which 
merely approximates the data collected from one series of experi- 

i 

ments or simulations if that function will provide a good fit to 
data from similar experiments or simulations. 
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The one-dimensional case will be considered here: Given 

a set of ordered pairs {(x^, y^)}, determine a function y - f(x) 
that corresponds to these ordered pairs in some fashion which can 
be called a ‘best fit.' The description is intentionally vague, 
at this point. Similar results can be obtained for n-dimensional 

spaces. The interested reader should consult intermediate to 

» 

advanced numerical analysis or approximation theory texts (example: 
Davis, 1965). The n-dimensional case can be proved using the 
techniques of Hilbert space. 

There exists a polynomial y ■ p(x) of degree at most n, which 

passes through any n + 1 points (x , y ), (x, , y,) ...» (x , y ) 

o o lx n n 

where x A < x. < . . . < x . It is not clear that p(x) is the best 
0 1 n r 

polynomial representation of the data. Consider Figure 1.1 
(Figure 4.8 from Conte and De Boor (1972)). It is exceedingly 
unlikely that any investigator would use a tenth-degree polynomial 
to fit this data, even though there is a tenth degree polynomial 
that interpolates it perfectly. Instead, a linear function was 
selected in Conte and de Boor. The investigator would assume that 
the small, seemingly random deviations from the straight line in 
Figure 1.1 represent error; sl*» would discount the possibility that 
the universe incorporates tenth-degree corrections to linear 
relationships. 

Linear regression and least-squares curve fitting are formal 
techniques for defining and obtaining a ’best fit' to data. Both 
methods attempt to fit data by minimizing the sum of the squares 




Least-squares straight-line approximation to 
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of the ’error,’ or distance from a data point to the best-fit curve. 
Multi-variable linear regression approximates the data via a linear 
equation in many unknowns. Curve-fitting approximates the data via 
a function, generally nonlinear, selected from a family of simple 
functions. A commonly-used family of functions is the polynomials; 
trigonometric functions or exponentials are sometimes used. Since 
the polynomials form a linear space, curve-fitting by polynomials 
is mathematically equivalent to multi-variable linear regression. 

One hardly needs mathematical formalism to deal with Figure 1.1. 
Consider, however, example 4.21 from Conte and de Boor (1972), 
presented in Table 1.1. Twenty-one data points are available. 

These data are to be fit by a polynomial. A 20th-degree polynomial 
can, in principle, be found to interpolate the data. The interpolating 
polynomials once obtained, would have a least-squares error of zero. 
However, the 20th-degree polynomial could not be found using simple 
algorithms on a computer. Computer truncation error would contaminate 
the result. Even if the 20th-degree polynomial could be obtained, 
it would be a poor choice as an approximant. The data in Table 1.1 

X 

happens to be. the values of e , rounded to the nearest 0.01. If 
this data were the result of a numerical simulation and the approxi- 
mating polynomial were to be used as a mathematical model of this 
growth, the approximant would be expected to resemble e x as closely 

X 

as possible, e is a very smooth function. The twentieth degree 

x 

interpolant to the rounded values of e is unlikely to be as smooth. 
Without computing the interpolating polynomial, it is not even 
possible to determine if it is monotonic on the interval (-1.0, 1.0). 



Table 1.1: Conte and De Boor (1972), Table 
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The interpolant would be a very poor choice of function for 
extrapolation outside the interval (-1, 1) . It would probably not 
represent the derivative of e X well. The twentieth degree inter- 
polating polynomial overfits the data; it faithfully reproduces 
the noise in the data so well that the signal tends to be obscured. 

What degree of polynomial best fits the data? The pragmatist 
would examine Figure 1.2, after Conti and ({>?. Boor (1972), Figure 4.9. 
Graphs (a) , (b) , and (c) are relatively smooth. Graph (a) certainly 
has a linear component, although it does exhibit some curvature. 
Graphs (b) and (c) resemble a quadratic and a cubic, respectively. 

In each case, the error graph resembles a polynomial of degree one 
higher than the degree of least-squares approximating polynomial 
whose error is being graphed. Graph (d) displays considerable 
oscillation, but the relative maxima at -1, 0.1, and 1 suggest a 
quartic component. Graphs (e) and (f) show virtually no regular 
pattern; they are pure noise. The fourth-degree polynomial (error 
graph (e)) fits the data as well as it can be fit. Note that 21 
data points were used to obtain a fourth-degree polynomial. Five 
points uniquely determine a quartic. 

The threshold at which overfitting begins depends on: 

1) the quality of the data 

2) The underlying function or physical law governing the 
data 

The more error the data contains, thr more likely that over- 
fitting will occur. To minimize the squares. the distance between 
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ORIGINAL PAGE *9 

OF POOR QUALITY 



The err<?r in the least-squares approximation to the data of Example 4.21 by polynomials 
of degree (a) zero, (b) one, (c) two, (d) three, (e) four, (f) five. 
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the data points and che approximating curve, the curve will be 
adjusted toward the data points. 1C the data points are scattered 
due to random errors, the curve may tend to oscillate away from the 
underlying function in order to approach the error-laden data points. 
By using a lower-degree polynomial as the least-squares approximant , 
this oscillatory behavior is suppressed. The number of local maxima, 
local minima, and inflection points of a polynomial is limited by 
the degree of the polynomial. The data is fit less well in this 
case, but the underlying function and its derivatives are fit 
better . 

Sensitivity studies on the original (Carlson and Boland, 1978) 
model indicate that M and P vary smoothly as functions of day and 
night temperature. Model output probably contains errors on the 
order of 1 K. Overfitting the model output by a polynomial is a 
real possiblility . Regression equations of high degree or regres- 
sions using many temperatures as predictors should be avoided. 
Another possible course of action, piecewise-linear interpolation 
in two variables, is presented in Appendix XI. 
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APPENDIX II: PIECEWISE LINEAR INTERPOLATION IN TWO DIMENSIONS 

II. 1 THEORETICAL OVERVIEW 

Piecewise-polynomial interpolation is an alternative to global 
least-squares regression for mathematically inverting the numerical 
model in the GOES method. The simplest case of piecewise- 
polynomial interpolation is piecewise-linear interpolation. 
Piecewise-linear interpolation has the advantage of simplicity 
over higher-order methods. However, to obtain an equivalent degree 
of fit to the data, additional data points will be needed. Thus, 
more than 16 (M, P) pairs may be needed in a piecewise-linear 
interpolation scheme. Model computation time is proportional to 
the number of (M, P) pairs. 

Consider a piecewise-linear interpolation scheme for diagnosing 
M from the two variables (predictors) DTD and DTIi. M can be con- 
sidered to be a function of the two variables DTD and DTN. If 
M(DTD, DTN) is multi-valued, this scheme and the regression scheme 
will fail. Therefore, assume that M(DTD, DTN) is single-valued. 
M(DTD, DTN) is a surface in (DTD, DTN, M) space. A piecewise- 
linear approximation to this surface consists of a number of 
planar segments lying "close” to the surface. A piecewise-linear 
interpolant to the surface at n specified (DTD, DTN) pairs is a 



105 


set of planar segments intersecting the M(DTD, DTN) surface at each 
of the n specified (DTD, DTN) values. If M(DTD, DTN) is not multi- 
valued, each planar segment can be expressed as an equation of the 
form: 

M = a^TD -I- b ± DTN + c (II. 1) 

The ith equation will be valid over the ith planar segment. 

Piecewise-linear interpolation in two dimensions, unlike 
global polynomial interpolation in one dimension, is not unique. 

Given a set of n (DTD, DTN, M) ordered triplets, where n > 3, 
there are, in general, several sets of planar segments interpolating 
the set of points.^ - A reasonable set of linear functions must be 
selected. Consider Figure II-l. Assume that the 16 points graphed 

in the (DTD, DTN) plane represented 16 data points from a model 

2 

run. Each point is associated with an M value, which is not plotted. 

I now address the problem of selecting the linear functions. 

Note that three points determine a plane. To find a piecewise- 
linear interpolant to the function M(DTD, DTN) over the domain 
defined by the points shown in Figure II. 1A, it is sufficient to 


Global polynomial interpolation in higher-dimensional spaces is 
unique, if the correct number of points is chosen. Linear inter- 
polation with two predictors is unique if three data points are 
chosen. 


2 

Sixteen data points from an actual model run are presented in 
Figure II-2. Note that the data points are close together at low 
values of DTD and DTN. The M and P values used in the model run 
were M and P values appropriate for the regression routine. To 
obtain a more uniform distribution of data points in the (DTD, DTN) 
plane for a piecewise-linear interpolation routine, different M 
and P values should be selected. 
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Figure II. 1: An example of triangulation. 
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Figure II. 2: Sixteen (DTD, DTN) pairs from a model simulation of 

the Kansas case (Section 4.3). Model initialization 
time was 0400 local time. Simulated GOES times were 
0800, 1400, and 2300 focal time. 
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break up that domain Into triangles and find planar functions 
over each triangle. Two ways of triangulating the region are 
presented in Figures II. IB and II. 1C. The pattern shown in 
Figure II. 1C is preferred, since the largest angle in each triangle 
is as small as possible. Once the region has been triangulated, 
the functions are relatively easy to find. At each vertex, 

(DTD, DTN, M) are known. Combine the three ordered triplets to 
obtain an equation of the form of Eq. II. 1. 

When this is done for all the triangles in the region, a 
piecewise-linear interpolant for M(DTD, DTN) over the domain will 
have been generated. Note that the interpolant is continuous at 
all points in the domain, but its derivatives fail to exist at the 
intersection of triangles. Since we have little interest in the 
derivatives of M with respect to DTD or DTN, we can accept the 
result . 

II. 2 PRACTICAL CONSIDERATIONS 

In the example <?f II. 1, the domain contained 16 (DTD, DTN) 
points. These points determined 18 triangles. The number of 
triangles needed to triangulate the domain increases faster than 
the number of data points included in the domain. We will assume 
that a fast algorithm for obtaining the coefficients a^» bj^ and ^ 
in Eq. II. 1 is available. The main practical problems are: 

1) triangulating the region efficiently 

2) given a (DTD, DTN) pair from satellite data, selecting 
the triangle containing the (DTD, DTN) pair quickly 
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The second problem cannot be finessed via lookup tables as was done 
In the GOES method and CD. Lookup tables could be used In those 
methods because: 

1) a global regression equation was used 

2) the variables DTD and DTN never appeared together In a 
terra In the global equation 

The piecewise-linear method meets the second criterion, but fails 
to meet the first. Since the domain is divided into irregular 
triangles and not rectangles, there is no hope of developing a 
simple lookup-table algorithm. 

The region could ba triangulated interactively, if desired. 
After the model has completed execution, the (DTD, DTN) data points 
could be displayed graphically. The operator could then enter 
triplets of points, defining triangles. The resulting triangles 
could then be displayed on the graphics terminal for verification. 
An option to delete points from the domain should be provided, 
since the model may generate several nearly coincident (DTD, DTN) 
pairs at high values of M and P. Another alternative would be 
automated selection of triangles. The automated selection could be 
preceded by an interactive opportunity to delete points from the 
domain, and followed by an interactive opportunity to revise the 
computer's choices. If a sufficiently robust triangulation routine 
can be developed, the entire triangulation procedure could be 


automated . 
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The efficiency of selecting the proper triangle (and thus the 
proper coefficients a^, b . , and c^) will depend on the data structure 
used to store the triangles, and the amount of storage space that 
can be used for redundant information. Since 16384 (DTD, DTN) 
pairs must be processed per image, storage should be sacrificed to 
improve efficiency. The triangle-selector and output-field builder 
should be combined: they should be in the same job, step. All 

desired output fields should be built concurrently, to avoid the 
necessity for searching more than once for the triangle containing 
a particular (DTD, DTN) value. For four output variables, 2K 
(2048 bytes) of storage will be needed to hold output fields. 

Storage requirements for the piecewlse-linear approximant should 
be modest: 12 coefficients plus four output variables, or 64 bytes 

per triangle. The data structure used to select the proper triangle 
can occupy most of the remaining core, if needed. 

The piecewise-1 inear inversion method should run considerably 
slower than the global regression-equation inversion method. No 
estimate is possible until the details of the triangle selection 
method are specified. Since the current inversion method takes on 
the order of 1-2 minutes per output field, a factor of five loss 
in efficiency is barely tolerable. There is reason to hope that 
the loss in efficiency will be less than a factor of five. 

Higher-degree piecewise-polynomial interpolants can be defined 
over a two-dimensional domain, but they have little or no 
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applicability to the problem at hand. Higher-order interpolants 
would produce better approximations to the output variables with 
fewer (DTD, DTN) pairs from the model. However, higher-order 
methods suffer from two flaws: 

1) One must be very careful to obtain a piecewise- 
polynomial approximation which is continuous at all 
points in the domain, and one which is never multi- 
valued 

2) The computation time required per (DTD, DTN) pair from 
the satellite data rises rapidly with increasing degree 
of polynomial approximant 

The first flaw probably can be overcome by careful programming. 

The second flaw appears •*<- be unremovable. Piecewise-polynomial 
interpolants of degree higher than one may be practical on main- 
frame computers or vector-processors, but they are not practical 
on minicomputers. 

II. 3 A PIECEWISE-LINEAR INTERPOLATION ALGORITHM 

The problem of selecting the triangle containing a particular 
(DTD, DTN) pair from satellite data is the key problem in the 
practical piecewise-linear algorithm. It may be possible to side- 
step this problem though the use of a very large lookup table, 
which would be built during interactive selection of triangles. 

The large amount of time needed to initialize the lookup table can 
be hidden from the user if asynchronous I/O or multitasking is 


used . 
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Under most circumstances, DTD and DTN will have ranges of 
less than 50 K. Since DTD and DTN have inherent uncertainties of 
+0.4 K due to the method of storing satellite Images on tape 
(Section 3.6), little additional uncertainty would be introduced 
into the method by mapping DTD and DTN into the integers from 1 
to 128. This mapping can be done by utility program MYC0NR2, 
already part of the GOES method library. Alternatively, a 
specialized version of KYC0NR2 could be written to map DTD and DTN 
into the integers with a minimum of additional error, using informa- 
tion about the original range of temperatures on the input (magnetic- 
tape) satellite images. The output from the specialized version 
of MYC0NR2 could be two-byte integers instead of one-byte Integers, 
to facilitate the use of the Integer values as indices into an 
array. The integer-DTD and integer-DTN fields would be written 
to disk as DEFINE FILE files, using the same record-numbering as 
is used to store the DTD and DTN fields (the record sizes would 
differ) . 

(DTD V , DTN) pairs generated by the numerical model would be 
displayed on the graphics terminal. Each (DTD, DTN) point would 
be represented by a distinguishing character and would be shown on 
the screen in its proper position in the DTD-DTN plane. While this 
display is being constructed, an auxiliary task would be clearing 
(setting to zero) a 128 by 128 array (16 K bytes) for use as the 
triangle-finding lookup table. This array will represent the 
region of DTD-DTN space of relevance to the GOES method. 
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After the array is cleared and all data points have been dis- 
played, the interactive selection of triangles would begin. At 
this point, asynchronous I/O or multitasking becomes essential. 

The following description assumes multitasking. Task A, the main 
task, will interact with the user and accept the triangvdation of 
the region. Task B, the subtask, will incorporate the selected 
triangulation into the lookup table. ^ 

Task A requests the user to input the first triangle by typing 
the distinguishing characters of the three vertices. Task A then 
plots the triangle on the graphics terminal, and asks the user to 
accept or reject it. When the first triangle is accepted, the 
multitasking begins. Task A attaches Task. B, and passes Task B 
the coordinates of the vertices of the first triangle. Task B, 
being able to count, knows that the vertices it has received are 
the vertices of the first triangle. Task B begins to fill all 
locations in the array corresponding to points on the edges and 
interior of this triangle with the integer 1 (one-byte integer, to 
save storage). When Task B has completed its work on triangle 1, it 
increments its triangle-counter, signals Task A, and goes to sleep. 

■^If asynchronous I/O were used, the actions performed by Task B would 
be interleaved with the terminal I/O in Task A. Task B actions 
would begin when the START I/O instruction to READ the terminal was 
issued. The terminal I/O event flags would be checked (and WAIT 
state entered, if necessary) on completion of Task B's work. Multi- 
tasking is preferred over asynchronous I/O, since multitasking 
simplifies the coding of efficiently-overlapping routines. 
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Meanwhile, Task A is asking the user to input the vertices of the 
second triangle. When the vertices are input, displayed, and 
accepted. Task A will: 

1) wait for Task B's signal, if it has not been received 

2) wake-up Task B and pass it the new set of vertices 

3) ask the user for the next triangle 

This procedure will continue until the region is completely 
triangulated. The user will inform Task A that the last triangle 
has been entered. Task A will then wait for Task B to finish 
setting the array locations representing the last triangle. When 
Task B signals. Task A will detach Task B. At this point, the 128 
by 128 array could be written to disk if it is to be caved. The 
array is now ready to be used to locate triangles. 

The linear equations for M, P, H q , and E can be computed by 
Task A or by Task B during the above procedure (let the faster 
routine do the extra work) , or they can be computed and stored 
after the triangulation has been completed (if storage is severaly 
limited). The index-numbers of the coefficients must match the 
index-numbers of the triangles used by Task B. 

The output fields can now be generated. Under the current 
image-storing convention, all. imagery used in the GOES method is 
stored in 128 by x.,3 arrays on disk. Each horizontal line of the 
image is stored in one logical record of a disk file; an image 
file consists of 128 data records and a control record. DTD, DTN, 
M, P, H q , and E are stored as REAL*4 variables: four bytes of 

storage per pixel. The integer-DTD and integer-DTN fields will be 
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INTEGER*2 fields: two bytes per pixel. All eight fields plus 


the lookup table must be 
storage requirements are 

simultaneously resident in core. 

The 

Lookup table: 

16 K 

16 K 

REAL*4 record: 

1/2 K each 

3 K 

INTEGER record: 

1/4 K each 

1/2 K 


Total storage for these fields: 19 1/2 K. 

(Note that only one record of each disk file must be resident at 
any one time, but that the entire lookup table must be resident at 
all times.) Additional storage for the piecevlse-linear coefficients 
will depend on the number of triangles used. Assume that 64 triangles 
are used — probably more than will ever be used in an actual case. 
For each triangle, assume that three coefficients are stored for 
each of the four output variables. These coefficients must be 
REAL*4 variables. The storage requirement for the coefficients is 
3*4*4*64 * 3072, or 3 K bytes. If the vertices of the triangles 
are to be stored, they will occupy under 1 K if stored as REAL*4, 
or under 1/2 K if stored as integers (coded as the integer-DTD 
and integer-DTN fields are coded) . This additional storage require- 
ment brings the total storage for all arrays to 23 K or 23 1/2 K 
of the vertices are kept; 22 1/2 K if they are discarded. The 
code required to perform the piecewise-linear approximation should 
easily fit in 4 K, allowing the job to run in a 28 K partition. 

(If the partition sizs is 32 K, an additional margin of safety 


exists) . 
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The output fields are produced by a simple algorithm. A 
(DTD, DTN) pair is selected from the temperature-difference input 
files. The corresponding integer-DTD, integer-DTN pair is also 
selected. The (integer-DTD, integer-DTN) pair is used as an index 
into the lookup table. The single byte returned from the lookup 
table is expanded into a two-byte integer via an EQUIVALENCE state- 
ment. If this integer is zero, the (DTD, DTN) pair represents 
out-of-range values and is treated accordingly. If the integer 
is nonzero, it is the sequence number of the triangle containing 
the (DTD, DTN) pair. The proper coefficients are obtained for each 
output variable from the tables of coefficients. The output values 
for the four output variables are obtained by substituting the 
(DTD, DTN) values into the appropriate equations of the form of 
Eq. II-l. 

Considerable care in coding will be necessary if this peicewise- 
linear interpolation and evaluation routine is to be implemented 
on a minicomputer with 32 K of storage available. The total 
storage requirement for the entire procedure is likely to exceed 
32 K. Assume that the routine is to be Implemented on the FSU 
Meteorology Department PDP 11/34 system. During the triangulation 
phase, the Grinnell driver and graphics routines will have to be 
resident. The DTD, DTN, M, P, H q , and E fields will not be needed. 
During the output-field generation phase, the various fields will 
have to be resident, but the graphics package will not be needed. 

The lookup table and the coefficient tables will be needed at all 
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times. An overlay structure will have to be devised to meet 
these requirements; appropriate sections of code must also be over- 
laid. It may also be necessary to force the system to allocate the 
absolute minimum input and output buffer space. In particular, 
it may be necessary to output the four computed fields without 
buffering, even at the cost of extra I/O time. Unbuffered I/O 
probably will not be needed if the available region size is 32 K 
(no system overhead) . 

If the piecewise-linear algorithm is implemented as described 
above, the execution time will not seem to be much longer than the 
execution time for the current regression algorithm. The evaluation 
step (done for each output variable) entails three table lookups 
(to get the coefficients), two multiply instructions, and two 
additions. One two-dimensional table-lookup is needed (to find the 
triangle) for each pixel; one fourth of this overhead can be assigned 
to the computation of each of the four output fields. The regression 
routine presently uses two table lookups and one addition per pixel 
per output field. The computation time for the piecewise-linear 
routine would be approximately three times that of the regression 
routine. There would, however, be a savings of 512 disk reads. 

The DTD and DTN fields will have to be read only once each, not 
four times (saving 768 reads) , and the new integer fields will 
have to be read (costing 256 reads). The initialization of the 
triangulation lookup table will be lengthy. However, the subjective 
computation time will be negligible. Most of the computation time 
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required for the initialization of the table will be concealed in the 
time spent in the interactive inputting of the triangles. The 
initialization will be done while the computer is waiting for the 
user to respond to requests for vertices. The user, busy thinking 
about vertices and typing on the keyboard, will be unaware of the 
time required for initialization of the large array. An automated 
procedure could be devised to select the triangles; however, if 
this were done, the user would become, painfully aware of the length 
of time required to initialize the large array. The user would also 
lose the option to delete points and the ability to select 
triangles of nonoptimal shape that are desirable for other reasons. 



